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Abstract 

Background: Electron Paramagnetic Resonance (EPR) is a non-destructive, non-invasive technique useful for the 
characterization of organic moieties in primitive carbonaceous matter related to the origin of life. The classical EPR 
parameters are the peak-to-peak amplitude, the linewidth and the g factor; however, such parameters turn out not to 
suffice to fully determine a single EPR line. 

Results: In this paper, we give the definition and practical implementation of a new EPR parameter based on the 
signal shape that we call the /?io factor. This parameter was originally defined in the case of a single symmetric EPR 
line and used as a new datation method for organic matter in the field of exobiology. 

Conclusion: Combined to classical EPR parameters, the proposed shape parameter provides a full description of an 
EPR spectrum and opens the way to novel applications like datation. Such a parameter is a powerful tool for future EPR 
studies, not only of carbonaceous matter, but also of any substance which spectrum exhibits a single symmetric line. 

Reproducibility: The paper is a literate program — written using Noweb within the Org-mode as provided by the 
Emacs editor — and it also describes the full data analysis pipeline that computes the /?io on a real EPR spectrum. 

Keywords: Electron paramagnetic resonance, Lineshape, Solid state chemistry, Carbonaceous matter, Exobiology, 
Literate programming, Python 



Background: Necessity for a shape factor definition 

In the field of exobiology, we need to determine the age 
of organic material in rock samples. Isotopic methods are 
commonly used to date the rock itself, but the organic 
matter may not be syngenetic with the rock. A novel solu- 
tion based on Electron Paramagnetic Resonance (EPR) 
was proposed [1]; it requires the determination of a new 
EPR parameter, the Rio, from the EPR spectrum of the 
rock sample, from which the age can be computed from an 
empirical log-linear correlation that was uncovered in [1]. 
Knowing the distribution of the different parameters that 
contribute to the Rio, we may also provide a confidence 
interval for the age thus determined. In the following, we 
shall explain what the classical EPR parameters are and 
what the proposed new parameter brings to the table, 
and then describe the algorithm for the determination of 
the Rio: how to process the data files generated during 
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an EPR experiment, extract the classical EPR parameters 
and compute their distribution in order to have an esti- 
mate of their error; compute the new Rio parameter and 
its distribution from the preceding distributions. Thanks 
to this paper, scientists may themselves extract the Rio 
parameter from EPR data and use it not only for datation 
purposes but also to uniquely characterize the observed 
EPR spectrum lineshapes. Our goal is to automate a man- 
ual process that has proved scientifically successful yet 
cumbersome and tedious when applied on datasets that 
are getting larger. In that version of our code, some of our 
algorithmic choices just mirror the — successful — manual 
process. We have chosen the Python language because of 
its high level, ease of development and popularity; last but 
not least, it also provides powerful libraries for scientific 
development, and speed of execution turned out not to be 
a key factor for our goals a . The Python code runs inside 
the Sage computing platform [2], which aims at providing 
a single computing environment both for numerical and 
symbolic computations. 
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Electron Paramagnetic Resonance (EPR) is a non- 
destructive and non-invasive technique which has indeed 
long been used for the study of paramagnetic defects 
(organic radicals) in carbonaceous materials. Such defects 
have been detected with high sensitivity in coals by pio- 
neering EPR works [3]. These types of radicals were 
therefore used for the characterization of a wide range of 
carbonaceous objects, ranging from coals [4-6] to cherts 
[7] through meteorites [8-11]. The EPR signal of kerogen 
is a single line, due to the presence of aromatic radi- 
cal moieties, with an unpaired electron spin delocalized 
in carbon p-type molecular orbitals [4,9,12,13]. Several 
parameters can be deduced from an EPR spectrum, based 
on the amplitude A pp , the line width AB pp and the res- 
onance field B res of the signal. However, for a single set 
of those three parameters, various lineshapes are possible 
(Figure 1); therefore, to fully determine the EPR line, a new 
EPR parameter, based on the lineshape, had to be defined. 

The shape of the magnetic resonance absorption line 
of a system of interacting and randomly distributed spins 
depends on the nature of the interactions (dipole-dipole 
or exchange), on the spin concentration and on the dimen- 
sionality of the spatial distribution of the spins [14-18]. 
This study is restricted to the case of a dipole-dipole 
type interaction between electron spins, thus excluding 
exchange interaction occurring in very concentrated elec- 
tron spin systems. Several limiting cases are distinguished 
in the literature, depending on the spin concentration and 
on the dimensionality of the distribution, cf. Table 1. 

In the high concentration regime (generally considered 
when the fractional site occupation r by a paramag- 
netic centre exceeds 0.1), the lineshape is approximately 
Gaussian [17]. This regime also occurs when the line 
is broadened by unresolved hyperfine interaction. Given 



Table 1 EPR lineshapes and lineshape parameter /? 10 for 
different limit regimes of dipolar broadening 
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Low: Lorentzian 


0 








2D 


Stretched Lorentzian 
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that EPR experimental spectra correspond to absorption 
derivatives, the Gaussian EPR line is described by: 
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where B is the applied magnetic field, B res the field at the 
centre of the line (maximum of absorption), A pp the peak- 
to-peak amplitude and Ai>pp the peak-to-peak linewidth 
(Figure 1). 

In the low concentration regime (generally considered 
when r < 0.01) with no hyperfine broadening, the 
lineshape depends on the dimensionality of the spatial 
distribution of the paramagnetic centres [16]. When the 
distribution is random, the resonance line may be calcu- 
lated from the relaxation function: 



G(d,t) =Bexp(-a-t^ 



(2) 



This function describes the decay with time t of the 
spin magnetization, perpendicular to the magnetic field, 
after an infinitely short microwave pulse. Parameter a is a 





^PP 


Gaussian 
d= 3 


/ ' ' 




d=2 

d- 1 


A5 PP 


' \ 




-* 1 

-4 -2 

(B- 


0 2 4 
-B res )/A5pp 


Figure 1 Theoretical EPR lines corresponding to upper limit cases of dipolar broadening. Continuous line: high spin concentration regime 
(Gaussian lineshape); Mixed line: diluted spin regime and 3D distribution (Lorentzian lineshape); Dashed line: diluted regime, 2D distribution 
(stretched Lorentzian); Dotted line: diluted regime, 1 D distribution (stretched Lorentzian). 
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constant that depends linearly on the spin concentration 
and parameter d represents the dimensionality of the spin 
distribution: d = 1 for a linear distribution, d = 2 for a 
distribution in a plane and d = 3 for a distribution in a 
volume. The EPR absorption is the Fourier transform of 
the relaxation function, and thus the EPR spectrum is the 
field derivative of this Fourier transform: 



F d (B-B res ) = m 

x exp 



-i(B — B res )t- 



h 



dt 



(3) 



where 9^ stands for the real part. In the case of a three 
dimensional distribution (d = 3), the EPR lineshape func- 
tion can be analytically calculated and corresponds to the 
field derivative of a Lorentzian function: 
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For lower dimension of spin spatial distribution (d < 3), 
the Fourier transform can only be calculated numerically. 
Figure 1 shows the theoretical EPR spectra corresponding 
to the Gaussian, Lorentzian (d = 3) and low dimensional 
(d — 1 and 2) cases. The wings of a Gaussian line fall 
off faster than those of a Lorentzian line while the wings 
of an EPR spectrum corresponding to a low-dimensional 
distribution fall off more slowly, giving rise to a so-called 
stretched Lorentzian lineshape. Originally, the R\§ line- 
shape factor was imagined after studying the spectra in a 
coordinate system (x,y) in which the difference between 
the lineshapes stands out more clearly [14], and where the 
Lorentzian becomes a straight line: 



JL(*)=*+ 7 ® 



and the Gaussian shape by an increasing exponential: 
/gW = exp(*- \) (6) 

with /g(x) > /lW; V#, cf. Figure 2. That coordinate 
system can be obtained thanks to the following transfor- 
mations as given in [14]: 
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where F = Fq or Fj. We shall thus define two functions, 
one that creates the new abscissas from the old x = B and 
the other that creates the new ordinates from the old x 
and y = F(B — B res ): 



(Functions to transform coordinates 3) = 
def xTransf orm (x, Bres, DeltaBpp) : 
return ( (x - Bres) / DeltaBpp ) 



(Ha) 



def yTransf orm (x, y, App, Bres, DeltaBpp) : 
return sqrt ( App * abs (x - Bres) / 
(DeltaBpp * abs (y) ) ) 

Following the Noweb literate programming style as 
described in [19], the above code is called a code chunk, 
with a unique name given between angle brackets and fol- 
lowed with an equal sign, together with a corresponding 
unique number made up of the page number and a letter 
starting at a and increasing alphabetically on a given page; 




■ Gaussian 
■d = 3 
■d = 2 
d=l 



Figure 2 Representation of the EPR spectra in the new (tfBenojBenc) coordinates system described by [14] and given in equation (7). 

Continuous line: Gaussian; mixed line: 3D distribution (Lorentzian); dashed line: 2D distribution (stretched Lorentzian); dotted line: 1 D distribution 
(stretched Lorentzian). 
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that number is mirrored in the left margin for easy refer- 
ence. The number on the end of line after the code chunk 
name indicates the code chunk where the current code 
chunk is used. Often, we shall add some code to an already 
existing code chunk, and that will appear in two different 
ways: first, the name between angle brackets will be fol- 
lowed by an equal sign attached to a plus sign (instead of 
a lone equal sign), and the numbers on the end of line will 
also indicate where the code chunk gets some new code 
(a small triangle is added to that number, i.e. < for previous 
existing definition, and > for the next new code). 

For diluted spin systems with low-dimensional dis- 
tribution, the representative function / lies below the 
line corresponding to a Lorentzian shape. To quantita- 
tively characterize the lineshape for systems intermediate 
between the above four ideal cases [Gaussian, Lorentzian 
(d = 3), one-dimensional (d = 1) and two-dimensional 
(d = 2)], we define a lineshape parameter measuring the 
deviation from a Lorentzian line as described in [7]: 

px=10 

R ™ = ^r\ [f(x)-f L (x)]dx (8) 

1U Jx=0 

This parameter corresponds to the algebraic surface 
between the curve / representing an experimental EPR 
spectrum and the curved representing a Lorentzian line. 
Rio is negative for a low-dimensional distribution (d < 3) 
and positive for an EPR line intermediate between 
Lorentzian and Gaussian lines (Table 1). The integration 
in equation (8) must be restricted to a finite range of x- 
values for the integral may not converge when x — »> 00. In 
practice, the range is limited to x < 10, since in most cases 
encountered the signal-to-noise ratio of the EPR spectra 
is poor for x > 10, inducing strong fluctuations in / and 
consequently in the lineshape parameter. Also, because 
of spectra with left/right assymmetry, the final R\o is the 
average of the values computed on the left and right of the 
resonance field, i.e. 

^10 = ^(^10+^-10) (9) 

To compute the integral in equation (8), we shall fol- 
low the method originally used: a simple top-left corner 
rectangular approximation. That allows full reproducibil- 
ity with the original manual method that was used before 
automation with a program; in the future we may replace 
it with a more accurate algorithm if there is a general 
agreement on the need to depart from the manual pro- 
cessing. We shall thus consider a matrix matrixXYL 
— a numpy array — made up of the abscissas of the 
spectrum in the first column, the ordinates of the spec- 
trum in the second column, and the ordinates of the 



ideal Lorentzian in the third column, with the number of 
lines corresponding to the number of data-points on the 
curves: 

4a (RIO rectangular integral 4a) = (5) 
Xspacing = (matrixXYL [:, 0] [1:] 

- matrixXYL [ : , 0] [ : -1] ) 
Ydifference = (matrixXYL [:, 1] [:- 1] 

- matrixXYL [ : , 2] [ : -1] ) 
Xspan = (matrixXYL [:, 0] [-1] 

- matrixXYL [ : , 0] [0] ) 

R10 = sum (Xspacing*Ydif f erence) 
/ Xspan 

The matrixXYL will be defined as a numpy array, and 
we use the sum function from the same library: 

4b (Import useful pylab functions 4b) = (12) 5b > 

from pylab import sum 

In order to construct the matrix matrixXYL, we 
need the data abscissas and ordinates and we use 
equation 5 for the yL coordinates of the ideal Lorentzian 
curve: 

5a (Matrix of XYL values in new coordinates 5a) = (5) 
matrixXYL = array ([ [xm, ym, yL] for x,y 

in zip (abscissas , ordinates) \ 
for xm,ym in [ (xTransf orm (x, Bres, 

DeltaBpp) , yTransf orm (x, y, App, 

Bres, DeltaBpp) ) ] 
for yL in [xm + 3/4.] if xm < 10 and 

(Same side of Bres 5c )] ) 

Again, we need to use the array data-structure, so we 
import it: 

5b (Import useful pylab functions 4b) + = (12) <4b 6f> 
from pylab import array 

Operationally, the R\q was only defined separately for 
the parts of the curve which abscissas x are larger or 
smaller than the resonance field Bres, and we thus define 
an operator testSameSideof Bres that will enable us 
to build two matrices matrixXYL, one for each side: 

5c (Same side of Bres 5c) = (5a) 
testSameSideof Bres (x,Bres) 

In the case of the left hand side, we look for x lower than 
Bres, and the opposite for the right hand side: 

5d (Left matrix of XYL values 5d) = (5h) 6a > 

testLef tSideof Bres = lambda x, 

Bres: x < Bres 
testSameSideof Bres = testLef tSideof Bres 
(Matrix of XYL values in new coordinates 5a ) 
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5e {Right matrix of XYL values 5e) = (5i) 
testRightSideof Bres = lambda x, 

Bres: x > 
Bres testSameSideof Bres = 

testRightSideof Bres 
{Matrix of XYL values in new coordinates 5a) 

We shall thus obtain two values of R\q, one for each side 

of£ r es> 

5f {Compute RIO 5f) = (lia) 
{RIO on the left 5h) 
{RIO on the right 5i> 
{Average RIO 5g) 

and we shall then use their average as the final value for 
the spectrum under study, cf. equation (9): 

5g {Average RIO 5%) = (5f) 
RIO = (R10_left + R10_right)/2 

5h {RIO on the left 5h) = (5f) 
{Left matrix of XYL values 5d) 
{RIO rectangular integral 4a) 
{Save RIO left value 5j) 

5i {RIO on the right 5i) = (5f) 
{Right matrix of XYL values 5e) 
{RIO rectangular integral 4a ) 
{Save RIO right value 5k) 

5j {Save RIO left value 5j) = (5h) 
R10_left = RIO 

5k {Save RIO right value 5k) = (5i) 
R10_right = RIO 

We need to be careful with the order of the values in the 
matrix giving the coordinates in the new coordinate sys- 
tem defined in equation (7): if we start from small values of 
x in the original frame, then, for the left hand side of B res , 
values in the new frame will decrease, whereas values on 
the right hand side will increase. Thus, values on the left 
side must be reversed, whereas that will not be necessary 
for the right hand side. 



6a {Left matrix of XYL values 5d)+ = (5h) <5d 

{Reverse matrix 6b) = 

6b {Reverse matrix 6b) = (6a) 
matrixXYL = matrixXYL [ : : - 1] 



Methods 

All the relevant discussion about the experimental part 
of the work, that involves collecting EPR data on the 



rock samples, can be found in [1]. In the current paper, 
we focus on the specific data handling and process- 
ing in order to extract the R\q parameter from an 
EPR spectrum and estimate the associated error. All 
computations were made in the Sage computing envi- 
ronment [2], with imports from the Numeric Python 
library [20]. 

In the spirit of reproducible research [21], the paper is 
written in the literate programming style [22]: the code 
and its explanation 13 are intertwined in a single place, and 
a particular program is then used to extract either the 
source code for execution on a computer or the liter- 
ate paper for reading by humans. Literate programming 
tools exist, and we use Noweb [19] and Org-mode [23,24] 
within Emacs with Evil mode to enable vi commands. We 
also make use of the Sagetex package that comes with 
the Sage distribution, that allows Sage code to be exe- 
cuted when compiling the LaTeX source of the paper 0 , 
and we have a home-built script that manages to combine 
Org-mode with Sagetex together with a Noweb output. 
Figures are produced either with Sage and Sagetex, or with 
Asymptote: it allows us to program figures, and thus make 
them executable, and embeddable in the LaTeX source 
code. The code will be made available through the teams 
website d . 

Processing data from an EPR file 

Removing the background signal 

EPR spectra on which the R\q factor was to be measured 
were selected for their symmetric and well-defined single 
absorption derivative signal. As usual in EPR studies, the 
large scale background signal was subtracted with a third 
degree polynomial fitted on the smooth parts of the spec- 
trum where the signal variations are only due to noise, 
which in practice correspond to the first and last 10% data 
points in a typical spectrum. 

6c {Remove background 6c ) = (12) 
{Compute number of data points in tails 6d) 
{Define points in tails 6e) 
{Compute polynomial on tails 6g) 
{Subtract polynomial from spectrum 7a) 

6d {Compute number of data points in tails 6d ) = (6c) 
numPointsTails = ceil (len (abscissas) /10 . ) 

6e {Define points in tails 6e) = (6c) 
tailX = concatenate 

( (abscissas [ : numPointsTails] , 
abscissas [-numPointsTails : ] ) ) 
tailY = concatenate 

( (ordinates [ : numPointsTails] , 
ordinates [-numPointsTails : ] ) ) 
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6f (Import useful pylab functions 4b) + = (12) <5b 6h> 
from pylab import concatenate 

6g (Compute polynomial on tails 6g) = (6c) 
degreePoly = 3 

backPoly = polyf it ( tailX, tailY, 
degreePoly) 

6h (Import useful pylab functions 4b) + = (12) <6f 7b> 
from pylab import polyfit 

From now on, the spectrum will be understood as the 
baseline corrected raw spectrum. 

7a (Subtract polynomial from spectrum 7a) = (6c) 
ordinates -= polyval (backPoly , abscissas) 

7b (Import useful pylab functions 4b) + = (12) <6h 7f> 
from pylab import polyval 

Reading the data for the spectra 

EPR Spectra are given as . txt files, with a name made up 
of the following informations: 

• sample-name 

• temperature-of-acquisition 

• microwave-power 

• number-of-scans 

For example, gunf lint_ambient_2mW_lscan . txt 
corresponds to a sample named gunf lint, studied at 
ambient temperature with a microwave power of 2mW 
using 1 scan e . 

7c (Load data 7c) = (12) 
(Define DATA directory 7g) 
(Define filename 7d) 
(Extract abscissas and ordinates 7e) 

7d (Define filename 7d) = (7c) 
fileName = 

' MB_gunf lint_ambient_2mW_lscan . txt ' 

The first two lines must be skipped when loading data: 
they provide the EPR acquisition parameters and the file 
description. EPR text files comprise three columns, giv- 
ing respectively the point index (starting from one and 
running to the total number of points recorded), the 
datapoint abscissa — the magnetic field B — and the dat- 
apoint ordinate —the intensity in arbitrary units. To ease 
data manipulation we extract two lists, abscissas and 
ordinates. 

7e (Extract abscissas and ordinates 7e) = (7c) 
data = loadtxt (DATA+ fileName , skiprows=2) 
#Skip first two comment lines 



abscissas = data [:, 1] . copy ( ) 
ordinates = data [ : , 2] . copy ( ) 

and the load function loadtxt will be taken from the 
pylab library. 

7f (Import useful pylab functions 4b) + = (12) <7b 8h> 
from pylab import loadtxt 

We also have to make sure that the DATA variable is 
defined, which is normally automatic within Sage: 

7g (Define DATA directory 7g) = (7c) 
try : 

DATA 

except NameError: 
DATA = 'data/' 

In order to plot the spectrum as in Figure 3, we use Sage 
builtin plot function list_plot. 

7h (Plot spectrum 7h) = (log 12) 

spectrumPlot = list_plot ( zip (abscissas , 

ordinates) ) 
spectrumPlot . save 

( DATA+ ' spectrumPlot .png' ) 

The distribution of the classical EPR parameters 

To uncover the underlying Lorentzian curve which will 
be compared to the original spectrum for the R\q com- 
putation, we need to find the three parameters that 
determine the latter: the peak-to-peak amplitude A pp , 
the line width AB pp and the resonance field i> r es- We 
define the peaks (positive and negative) as the extrema 
of the spectrum ordinate values, and the A pp and AB pp 
as the difference between the peaks' ordinates and 
abscissas, respectively. 

8a (Compute App 8a) = (8d) 
Amax = max (ordinates) 
Amin = min (ordinates ) 
App = Amax - Amin 

8b (Compute DeltaBpp 8b) = (8d) 
indexAmax = list (ordinates) . index (Amax) 
Bmin = abscissas [indexAmax] 
indexAmin = list (ordinates) . index (Amin) 
Bmax = abscissas [indexAmin] 
DeltaBpp = Bmax-Bmin 

The resonance field Bres was defined as the value at 
which the EPR lineshape crosses the baseline of the spec- 
trum, which corresponds to the zero axis since the spectra 
are baseline corrected. 
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TO 



Figure 3 The loaded EPR spectrum (dots) and the corresponding theoretical Lorentzian (continuous): the Rio factor is based on the 
integral difference between the two, cf. equation (8). 



8c {Compute Bres 8c) = (8d) 
for y in 

(The ordinates between the two extrema 8e): 
(Find when ordinate crosses baseline 8f) 
Bres = 

(Mean of the two ordinates above and below 
baseline 8g) 

8d (Compute the classical EPR parameters 8d) = (9g) 
(Compute App 8a) 
(Compute DeltaBpp 8b) 
(Compute Bres 8c) 

8e ( The ordinates between the two extrema 8e) = (8c) 
ordinates [indexAmax : indexAmin+1] 

8f (Find when ordinate crosses baseline 8f) = (8c) 
if y < 0 : 

yCross = y 
break 

The resonance field Bres is thus the mean of the two 
ordinates lying above and below the baseline respectively: 



8h (Import useful pylab functions 4b) + 
from pylab import mean 



(12)<7f9b> 



(Mean of the two ordinates above and below 
baseline 8g) = 
mean (map (lambda v: abscissas [list 
(ordinates) . index (yCross) + v] 
[-1, 0])) 



(8c) 



Knowing the distributions of the classical EPR param- 
eters App, DeltaBpp and Bres, we may check visually 
their normality thanks to a histogram plot; if normal, we 
may propagate their standard deviation in the global R\q 
error calculation. 

9a (Plot histograms of the classical EPR 

parameters 9a) = (9g) 
setEPRparams = [listApp, listDeltaBpp , 
listBres] 

for i in range (len (setEPRparams) ) : 
elf () 

hist (setEPRparams [i] ) 

savef ig (DATA+ ' distrib' + str(i)) 

9b (Import useful pylab functions 4b) + = (I2)<8h9d> 
from pylab import elf, hist, savef ig 

9c (Compute the moments of the classical EPR 

parameters 9c) = (9g) 
[App, DeltaBpp, Bres] = 

map (mean, setEPRparams) 
[sigmaApp, sigmaDeltaBpp , sigmaBres] = 
map (std, setEPRparams) 

9d (Import useful pylab functions 4b) + = (I2)<9b9f> 
from pylab import std 
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In order to uncover the classical EPR parameters' dis- 
tributions, we chose the Monte Carlo error propagation 
method, cf. [25]: we take the measured spectrum, consider 
each data point as the mean of a random variable, then 
draw a new value for each data point given its distribu- 
tion. For that, we suppose it is a normal distribution, with 
mean given by the data point and standard deviation given 
by the square root of the mean f ; we thus use the normal 
distribution generator provided by randn in the pylab 
library. 

9e {Add noise to data 9e) = (9g) 
ordinates += sqrt (abs (ordinates) ) 

* randn (len (ordinates) ) 

9f {Import useful pylab functions 4b) + = (12) <9d 

from pylab import randn 

With this approach, a large number of cloned data sets is 
generated, for which App, Del taBpp and Bres are com- 
puted; we then check for their normality by plotting their 
distribution and, if confirmed, compute their standard 
deviation for later use when computing the distribution of 
the £10. 

9g {Distribution of the classical EPR parameters 9g) = 

(12) 

{Create the lists for the classical EPR parameters 9h ) 
{Backup dataiOc ) 

{Repeat a large number of times 10b) 

{Add noise to data 9e) 

{Compute the classical EPR parameters 8d) 

{Store the parameters 9i) 

{Retrieve original data iod) 
{Plot histograms of the classical EPR parameters 9a) 
{Compute the moments of the classical EPR 

parameters 9c) 

To store the parameters, we need to create the three 
empty lists listApp, listDeltaBpp and listBres. 

9h {Create the lists for the classical EPR 

parameters 9h) = (9g) 
listApp, listDeltaBpp, listBres 
= [] , [] , [] 

We then use the append function to add each calcu- 
lated set of data to the storage lists. 

9i {Store the parameters 9i) = (9g) 
listApp . append (App) 
listDeltaBpp . append (DeltaBpp) 
listBres . append (Bres) 
/ 



For the Monte Carlo error propagation, we need to iter- 
ate a sufficient number of times in order to produce a 
significant set of data ; we thus create a global variable that 
specifies the number of Monte Carlo iterations. 

10a {Define global variables 10a) = (12) 
numlter = 100 

10b {Repeat a large number of times 10b) = (9g lia) 

for _ in xrange (numlter) : 

Because we add some noise during the Monte Carlo 
error propagation, and thus modify the original data, we 
need to store it before starting the Monte Carlo and 
retrieve it for each iteration in the Monte Carlo. 

10c {Backup data lOc) = (9g) 
originalOrdinates = ordinates . copy ( ) 

iod {Retrieve original data iod) = (9g) 
ordinates = originalOrdinates . copy ( ) 

Extracting the new Rio factor from the spectrum 

The R\q factor is calculated from the difference with the 
ideal Lorentzian derivative, which equation is: 

-A pp ABl p (x-B res ) 
(|A£2 p + (*-£ res ) 2 ) 2 

where A pp is the signal amplitude, AB pp the siglnal width 
and B res the resonance field; such an expression supposes 
that the background signal has been subtracted, i.e. that 
A moy = 0. We thus compute the theoretical Lorentzian 
ordinates yL corresponding to the same abscissa as that of 
the spectrum and the same classical EPR parameters App, 
DeltaBpp and Bres as that of the spectrum; we store 
them in a list lorentzOrdinates. 

i0e {Define theoretical Lorentzian i0e) = (lOf) 
lorentzFun = lambda x: -App*DeltaBpp^3 
* (x-Bres) / (0 . 75*DeltaBpp^2 
+ (x-Bres) x 2) ^2 
lorentzOrdinates = [lorentzFun (x) 
for x in abscissas] 

We plot the spectrum and its corresponding Lorentzian 
curve for visual checking. 

lOf {Plot Lorentzian lOf) = (log) 
{Define theoretical Lorentzian i0e) 
lorentzPlot = list_plot ( zip (abscissas , 

lorentzOrdinates) ,\ 
colors red', plot j oined=True) 
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lorentzPlot . save 

( DATA+ ' lorentzPlot .prig' ) 

log {Plot spectrum and Lorentzian log) = (12) 
{Plot spectrum 7h) 
{Plot Lorentzian lOf) 

dlPlot = spectrumPlot + lorentzPlot 
dlPlot . save (DATA+ ' dlPlot . png' ) 

Now the 7?io parameter is computed relatively to the 
theoretical Lorentzian having the same set of classical EPR 
parameters, so we could compute the error on the for- 
mer by propagating analytically the errors of the latter, 
which we now know thanks to the previous application 
of the Monte Carlo error propagation method. However, 
we found it easier and somewhat more in line with the 
computational approach to use a Monte Carlo approach 
to propagate the errors. We thus need to repeat the R\$ 
computation for a series of values of Bres, DeltaBpp 
and App to which we add a random error compatible with 
their distributions^ 

lia {RIO distribution computation lia) = (12) 
{Functions to transform coordinates 3) 
{Create RIO list lie) 
{Backup classical EPR parameters lib) 
{Repeat a large number of times 10b) 

{Add noise to classical EPR parameters lid) 

{Compute RIO 5f) 

{Store RIO llf) 

{Retrieve original classical EPR parameters lie) 
{Compute mean and standard deviation 
for RIO ii g ) 

Because we modify the classical parameters during the 
Riq computation, we need to store the values and retrieve 
them before and after each iteration of the Monte Carlo: 

lib {Backup classical EPR parameters lib) = (lia) 
backupEPR = [App, DeltaBpp, Bres] 

lie {Retrieve original classical EPR parameters lie) = 

(Ha) 

[App, DeltaBpp, Bres] = backupEPR 

lid {Add noise to classical EPR parameters lid) = (lia) 
App += randn ( ) *sigmaApp 
DeltaBpp += randn () *sigmaDe It aBpp 
Bres += randn () *sigmaBres 

lie {Create RIO list lie) = (lia) 
listRlO = [] 



iig {Compute mean and standard deviation 

for RIO 11%) = (lia) 
RIO = mean (listRlO) 
sigmaRlO = std(listRlO) 

Results and conclusion 

We now have extracted the R\$ parameter together with 
its distribution and may proceed to use it, for example to 
determine the age of organic matter inside rock samples 
[1]. Given the distribution, we may then check if the mean 
and standard error do indeed properly characterize the 
parameter, and eventually assign a probability to a range 
of ages for the rock sample. The code runs in only a few 
minutes, if we take into account all the Monte Carlo com- 
putations. In [1], we demonstrate that the data processing 
as reported here can indeed provide us with a reasonable 
estimate for the age of rock samples older than 1 billion 
years. 

The complete code 

12 {rlO.py 12) = 

{Import useful pylab functions 4b) 
{Define global variables 10a) 
{Load data 7c) 
{Remove background 6c) 
{Plot spectrum 7h) 

{Distribution of the classical EPR parameters 9g) 
{Plot spectrum and Lorentzian lOg) 
{RIO distribution computation 11a) 

Endnotes 

a Anyway, tools exist to go faster when needed, as 
Cython inside Sage that allows easy variable typing. 

b Or maybe the explanation and its code. . . literate 
programming is really a whole new approach to writing, 
thinking and coding. 

c This means that the outputs of some code need not be 
pasted inside the paper, but can be computed on the fly 
as needed. 

d The url is http://hpu4science.org. 

e This sample is part of the study where the R\o 
parameter was proposed as a datation method [1]. 

f This corresponds to a normal distribution arising from 
a Poisson distribution, and is the common practice in 
EPR because of the underlying counting process when 
measuring the absorption giving the spectrum. We can 
indeed check it is so by studying the noise on the flat tails 
of EPR spectra. 

g Using the Monte Carlo approach would also allow us 
to draw the values for the classical parameters according 
to their computed distribution. 



llf {Store RIO llf) = 

listRlO . append (RIO) 



(Ha) 
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Chunks 

(Add noise to classical EPR parameters 1 1 d) 
(Add noise to data 9e) 
(Average RIO 5g) 

(Backup classical EPR parameters 1 1 b) 
(Backup data 1 0c) 
(Compute App 8a) 
(Compute Bres 8c) 
(Compute DeltaBpp 8b) 

(Compute mean and standard deviation for RIO 1 1 g) 
(Compute number of data points in tails 6d) 
(Compute polynomial on tails 6g) 
(Compute RIO Sf) 

(Compute the classical EPR parameters 8d) 

(Compute the moments of the classical EPR parameters 9c) 

{Create RIO list 11e) 

(Create the lists for the classical EPR parameters 9h) 

(Define DATA directory 7g) 

(Define filename 76) 

(Define global variables 10a) 

(Define points in tails 6e) 

(Define theoretical Lozentzian 1 0e) 

(Distribution of the classical EPR parameters 9g) 

(Extract abscissas and ordinates 7e) 

(Find when ordinate crosses baseline 8f) 

(Functions to transform coordinates 3) 

(Import useful pylab functions 4b) 

(Left matrix of XYL values 56) 

(Load data 7c) 

(Matrix of XYL values in new coordinates 5a) 

(Mean of the two ordinates above and below baseline 8g) 

(Plot histograms of the classical EPR parameters 9a) 

(Plot Lorentzian 1 0f) 

(Plot spectrum 7h) 

(Plot spectrum and Lorentzian 10g) 

(RIO distribution computation 11a) 

(RIO on the left 5h) 

(RIO on the right 5i) 

{RIO rectangular integral 4a) 

(rW.pyU) 

(Remove background 6c) 

(Repeat a large number of times 10b) 

(Retrieve original classical EPR parameters 1 1 c) 

(Retrieve original data 10d) 

(Reverse matrix 6b) 

(Right matrix of XYL values 5e) 

(Same side of Bres 5c) 

(Save RIO left value 5j> 

(Save RIO right value 5k) 

(Store RIO 11f) 

(Store the parameters 9i) 

(Subtract polynomial from spectrum 7a) 

(The ordinates between the two extrema 8e) 
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